function [ resonantFreq ratio freq idx resp stim respf stimf ] = ZAP( prefix, traceNums, minFreq, maxFreq, channelNum)
%ZAP Calculates the impedance amplitude profile from the response to a
%chirp stimulus
%
% Calculates impedance as the magnitude of the ratio of the Fourier
% transform of the response (in mV) to the Fourier transform of the 
% stimulus (in pA). See Shin 2008, Puil 1986. The maxFreq arguments only
% applies to the figure and resonantFreq search (the transforms contain all
% frequencies sampled).
%
% Arguments:
%  prefix       path to cell, e.g. 'BS0065/BS0065AAAA'
%  traceNums    which traces to include, e.g. 42 or [ 42:49 ]
%  minFreq      min frequency (in Hz) to include in plot and search
%  maxFreq      max frequency (in Hz) to include in plot and search for 
%               resonant frequency, e.g. 20
% 
% 2009-09-28 Ben Suter
% 
% add another argument for channel 2. Channel number (1 or 2) has to be
% chosen.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TODO: Would it be better to calculate the impedance for each trace
% separately and average in the end?

% TODO: Conceptual problem: Frequencies are not all treated equally in the
% stimulus: some are positive injected current, some are at negative
% injected current (or, some cause more, some less polarization of the
% cell). The membrane is active, so can this periodic asymmetry may lead to
% periodic artifacts in the impedance profile? I see some ... perhaps the
% solution is to lower the current amplitude of the stimulus, to get about
% a +/- 5 mV voltage response (small around resting potential). I've made
% such a stimulus, but we may need finer control. Make a user function to
% set the amplitude, and write to autonotes, and make a hotswitch.
% Note: GS said to avoid causing the cell to "spike", but I think that is
% not stringent enough - because of the periodic asymmetry I describe
% above.
% If I'm right, then a lower amplitude stim should result in the same or
% similar impedance profile, but without the periodic artifact!

for f=1:length(traceNums)
    xsg = [ prefix sprintf('%04.0f', traceNums(f)) '.xsg' ];
    chirp = load(xsg, '-mat');
    sampleRate(f,:) = chirp.header.ephys.ephys.sampleRate;
    if channelNum == 1
        resp(f,:) = chirp.data.ephys.trace_1;
    else
        resp(f,:) = chirp.data.ephys.trace_2;
    end
    stimTemp = chirp.header.ephys.ephys.pulseParameters{1};
    stimTemp = stimTemp.signal'; % make same dims as response
    stim(f, :) = stimTemp(1:end-1); % has one data point more than the response
end
if numel(unique(sampleRate)) ~= 1
    error('All traces must use the same sample rate');
else
    sampleRate = sampleRate(1);
end
resp = mean(resp, 1);
stim = mean(stim, 1);

freq = (0:length(resp)-1) * sampleRate / length(resp);
respf = fft(resp);
stimf = fft(stim);
% Impedance = voltage / current
% response is in mV, stim is in pA, so ratio is in GOhm
ratio = abs(respf ./ stimf);

% Plot impedance vs. frequency
% Omit the first data point (0 Hz) and stop after maxFreq
% idx = 2: maxFreq / sampleRate * length(resp) + 1;
idx = round(minFreq / sampleRate * length(resp)) + 1 : round(maxFreq / sampleRate * length(resp)) + 1;
figure;
plot(freq(idx), ratio(idx) * 1000); % display in MOhm
title({'Impedance Amplitude Profile (ZAP) for', prefix});
xlabel('Frequency (Hz)');
ylabel('Impedance (MOhm)');

[mx ix] = max(ratio(idx));
resonantFreq = freq(idx(ix));

% Draw a vertical line indicating the resonant frequency
x = [resonantFreq resonantFreq];
y = ylim(); y(2)= mx*1000;
line(x, y, 'Color', 'r', 'LineStyle', ':', 'Marker', 'o');
 
end

